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Abstract 

I The time evolution of several interacting Ginzburg-Landau vortices according to an equa- 

' tion of Schrodinger type is approximated by motion on a finite-dimensional manifold. That 

Tlj- ■ manifold is defined as an unstable manifold of an auxiliary dynamical system, namely the 

[ gradient flow of the Ginzburg-Landau energy functional. For two vortices the relevant un- 

f-*) I stable manifold is constructed numerically and the induced dynamics is computed. The 

CN ■ resulting model provides a complete picture of the vortex motion for arbitrary vortex sep- 

, aration, including well-separated and nearly coincident vortices. 

« ■ 

I AMS classification scheme numbers: 35Q55, 37K05, 70K99 

> 

^ ! 1 Ginzburg-Landau vortices and their dynamics 

. 

Vortices play a fundamental role in a large variety of physical systems, ranging from ordinary 
fluids over condensed matter to the early universe. This variety is reflected in the mathematical 
models used to describe the formation, structure and dynamics of vortices. In fluid dynamics, 
the basic ingredient of the mathematical model is the velocity field of the fluid, and the vortex is 
a particular configuration of that velocity field. In condensed matter theory, the mathematical 
models are field theories. The basic field of such models is a complex valued scalar field, and 
vortices are particular configurations of that field. One important difference between vortices 
in ordinary fluids and those in condensed matter systems is that the total vorticity can take 
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arbitrary real values in the former but is quantised in the latter. Nonetheless, it is sometimes 
possible to establish a precise mathematical connection between field theory models and those 
used in ordinary fiuid dynamics. 

In this paper we are concerned with vortices in the Ginzburg-Landau model. For general 
background and references we refer the reader to the book pQ . We are interested in the dynamics 
of such vortices in the plane, so the basic field of the model is the function 

^/'iMxR^^C, (1.1) 

depending on time t G M and spatial coordinates x = {xi, X2) € M^. Sometimes we also use polar 
coordinates (r, 6) in the plane and parametrise the field ip in terms of its argument and modulus 
via = ^/pexp{i(j)), where p = l^/'p. In this paper ip will be used to denote both time-dependent 
and static configurations. For static configurations we often suppress the first argument and 
simply write iIj{x). For configurations varying with time we sometimes write ip for the partial 
derivative dtip. Often we write di for d/dxi, dt for d/dt and so on. 

The Ginzburg-Landau energy is the following functional of iJj: 

m = \j v^v^ + \m'~ If d'x. (1.2) 

We impose the boundary condition 

hm\^{r,e)\ = l, (1.3) 

r— >oo 

which implies that for sufficiently large R the field 

MO) = ^^iR,e) (1.4) 

is a well-defined map from the spatial circle of large radius R to the unit circle in C. It therefore 
has an associated integer winding number or degree, which can be computed via the following 
line integral over the circle S}^ of radius R 

deg[ip] = ^— ® dlnip. (1.5) 

This integer, denoted n in the following, is also called the vortex number. 

The variational derivative of the Ginzburg-Landau functional ()1.2|) with respect to ip is 

H = -1(A^ + (1-|^|2)^) (1.6) 
and, because of the reality of E, 

^^=H. (1,7) 
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The Ginzburg-Landau equation is the equation for critical points of the Ginzburg-Landau energy, 
so it reads 



-A^ + dV-p- 1)^ = 0. (1.8) 

A Ginzburg-Landau vortex is a solution of this equation with non-vanishing vortex number n. 
It is shown in j2] that for any configuration of non-vanishing vortex number which attains the 
limit (jl.Hj) uniformly in 9 the Ginzburg-Landau energy is necessarily infinite. The origin of this 
divergence is not difficult to understand and also explained in The essential point is that 
the gradient terms in the energy density contain the term \doil)\^ /r"^ which, for n 7^ 0, leads to a 
logarithmic divergence. Since the divergence only depends on n and not on any other details of 
the configuration, it can be removed by introducing a smooth cut-off function 



^^(^) = \ 1 io.\x\>R + R-^ ^^-9) 
and defining the renormalised energy functional 

E..A^] = \ I V^l^V^^-^^^XR + lm'-^fd'x. (LIO) 

This renormalisation procedure is natural both from the point of view of physics and in numer- 
ical investigations. The point is that the total vortex number n is conserved during the time 
evolution. When studying the interacting dynamics of n vortices during a finite time interval we 
can choose R so large that all the vortices remain well inside the disc of radius R during that 
time interval. The renormalisation procedure removes a divergence which only depends on the 
conserved quantity n but not on any other details of the dynamics. 

In studying the dynamics of vortices one has a choice of several different ways of extending the 
static Ginzburg-Landau equation to a time-dependent evolution equation. We are interested in 
the following first-order equation of Schrodinger type, often called the Gross-Pitaevski equation: 

^^ = -Ai; + m'-l)^. (Lll) 

The vortex dynamics dictated by this equation has been studied in a large number of publica- 
tions. In two seminal papers [SI 1^ Neu showed that in a scaling limit where the vortex size 
shrinks to zero, the time evolution according to the partial differential equation reduces to 

a set of coupled ordinary differential equation for the centres of vorticity. Moreover, he showed 
that set of ordinary equations to be the Kirchhoff-Onsager law for the motion of fluid vortices in 
incompressible, nonviscous two-dimensional flows. Neu's work has in turn inspired a number of 
authors. Some, mathematically motivated, have investigated the scaling limit further and have 
put his work on a more rigorous mathematical footing, see [312] and [7j. Others, starting with 
the interpretation of Neu's work as a flnite-dimensional approximation to the inflnite-dimensional 
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dynamical system governed by the Gross-Pitaevski equation, have tried to go beyond this approx- 
imation. Physically, one may think of Neu's approximation as a point-particle approximation to 
vortex dynamics in a field theory. This approximation is expected to be reasonable when the 
vortices are well-separated and moving slowly. However, when the vortices overlap, the point- 
particle approximation is poor. Similarly, when the vortices move rapidly we expect them to 
excite radiation in the field theory which is not captured by the point-particle approximation. 

More recently, Ovchinnikov and Sigal have derived Neu's approximation from a different point 
of view and have computed leading radiative corrections to it in a series of papers IHUHl- They 
work with the Lagrangian from which the Gross-Pitaevski equation can be derived and explicitly 
study the truncation of the infinite dimensional Gross-Pitaevski dynamical system to a finite- 
dimensional family of multi-vortex configurations with pinned centres of vorticity. The induced 
Lagrangian on that finite dimensional family reproduces the Kirchoff-Onsager law when the 
vortex centres are well-separated. Moreover, the approach makes it possible to study the coupling 
between the vortex motion and radiative modes and to compute radiative corrections. However, 
Ovchinnikov and Sigal's family of pinned multi-vortex configurations is less useful when trying to 
understand the dynamics of overlapping vortices. When vortices get close together pinning the 
vortex centres is mathematically awkward and physically unnatural. In this paper we study the 
truncation of the Gross-Pitaevski dynamical system to a finite dimensional dynamical system 
using a different set of multi-vortex configurations. Our configurations are similar to those of 
Ovchinnikov and Sigal when the vortices are well-separated, but we believe they provide a more 
accurate description of vortex dynamics when the vortices overlap. Our approach is inspired by 
an approximation scheme proposed by Manton in the context of Lagrangian soliton dynamics. 
Manton considered the problem of defining a smooth finite-dimensional family of multi-soliton 
configurations which could be used as the configuration space for a finite-dimensional low-energy 
approximation to multi-soliton dynamics. His proposal is to consider an auxiliary evolution 
equation, namely the gradient fiow in the potential energy functional, and to use the unstable 
manifold of a suitable saddle point as the truncated configuration space. Such an unstable 
manifold is the union of paths of steepest descent from the saddle point. In practice this scheme 
is not easy to implement, and so far it has only been used to study soliton dynamics in one spatial 
dimension JTT\. In this paper we show that it is very well suited to studying the dynamics of 
overlapping Ginzburg-Landau vortices. The basic reason is that in the Ginzburg-Landau model 
there are well-known rotationally symmetric saddle point solutions for all vortex numbers |n| > 1 
which can be used for the construction of the unstable manifold. We focus on the case n = 2 and 
find that the relevant unstable manifold is two-dimensional. We describe its geometry, compute 
the induced Lagrangian and use it to study the relative motion of two vortices at arbitrary 
separation. 
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2 The Gross-Pitaevski equation and its symmetries 



For our study of the Gross-Pitaevski equation it is essential that we can derive it from a La- 
grangian. The required Lagrangian is the following functional of time-dependent fields ip 

m=T[^,i;]-E[i;], (2.1) 

where the kinetic energy functional T is given by 

T{^, iP) = j lm{^piJ) d^x (2.2) 

and the potential energy functional E is the Ginzburg-Landau functional ()1.2j] . One checks that, 
essentially because of the linearity of the kinetic energy in ip, the total energy or Hamiltonian is 
equal to the potential energy E. Using the formula ()1.6|) . the Euler-Lagrange equation of L 

SE 

dttP = -2tj^, (2.3) 

is readily seen to be the Gross-Pitaevski equation 

We define the configuration space to be the space of smooth fields satisfying the boundary 
condition ()1.3|) . i.e. 

C = : M2 ^ C| lim |V(a;)| = 1}. (2.4) 

As mentioned earlier, the boundary condition means that every ip & C has an associated integer 
winding number or degree. Using Stokes's theorem and taking the radius R to infinity in p.5|) 
one shows that the degree can also be written as the integral 



degM = / 7, (2.5) 

where the integrand 

7 = —dtp A dtp (2.6) 
27ii 

is called the vorticity. The integer degree cannot change under the continuous time evolution, 
which can therefore be restricted to one of the topological sectors 

Cn = {i^e C\ deg[ij] = n}. (2.7) 

For each n, the pair (C„, L) is an infinite dimensional Lagrangian dynamical system. From the 
Lagrangian viewpoint C„ is the configuration space and a motion is a path t \—>- tp{t, ■) in Cn 
which satisfies the evolution equation Alternatively, we can adopt the Hamiltonian point 

of view. Then C„ should be thought of as the phase space and the kinetic part of the Lagrangian 
(|2.H) defines a symplectic structure on C„. This is a two-form Q on the tangent space of C„. For 
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a more careful discussion of that tangent space we refer the reader to [2], but for our purpose it is 
sufficient to think of tangent vectors as complex valued functions which tend to zero at infinity. 
If ^ is such a tangent vector we use the notation introduced in |2j for elements of the complexified 
tangent space: 

f = (I) ■ (2.8) 



A, 

If ^ and rj are two elements of the same tangent space the symplectic form is defined via 

^i^V) = ll j{in-lri)d\. (2.9) 
It follows that the Poisson bracket of ip and if) is 

{^,7/>} = 2z. (2.10) 
One checks that the Gross-Pitaevski equation can then be written in the canonical form 

^ = {E[i^ii^}. (2.11) 

As an aside we point out that the total vorticity or degree is related to the pull-back of the 
symplectic form VL evaluated on the vector fields di and 82 on : 



r^{d,,d2) = ~^ / 7- (2.12) 

Thus ii ip & Cn then ilj*VL{di, 82) = —mi. 

The Lagrangian L has a large invariance group. Writing R G 5*0(2) for spatial rotations and 
G for translations in the plane, elements {R, d) of the Euclidean group 

E2 = S0{2) tx R2 (2.13) 

act on fields via pull-back 

ip ^ ip o (R, d)'\ (2.14) 

i.e. ip o [R,d)~^{x) = iIj{R~^{x — d)). Both the Lagrangian and the degree are invariant under 
this action. Similarly, unit complex numbers G U{1) act on fields ip as phase rotations 

%Ij ^ e*"V (2.15) 

and leave both the Lagrangian and the degree invariant. Spatial refiections about any axis leave 
the Lagrangian invariant but change the sign of the degree. The same is true for the combination 
of complex conjugation with time reversal T : t ^ —t . However, with S : {xi,X2) ^ (xi, —X2), 
the combination 

C ■ ip ^ %l) o {ST) (2.16) 
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leaves the degree and the Lagrangian invariant. To sum up, for each n G Z the group 

C K (^2 X f/(l)) (2.17) 

acts on Cn preserving the Lagrangian L and hence the equations of motion. The two-dimensional 
Galilean group provides an additional more subtle symmetry. If iplt, x) solves the Gross-Pitaevski 
equation, then so does the configuration 

x) = e*(^"-''-^'*V(^, x-vt), (2.18) 

where the parameter w e is physically interpreted as the boost velocity. 

To end this section we note the conservation laws which follow from the symmetry group 
fl2.17|l . The Noether charge which is conserved due to the invariance under phase rotations is 

g[^]= / \^P\'£x. (2.19) 

Invariance under spatial rotations leads to the conserved charge 

J[iP] = I \m{ipdei)) d^x (2.20) 

and invariance under translations leads to the conservation of the vector (Pi, P2) with components 

Pi[^] = / lm{ipdi^lj)d^x i = l,2. (2.21) 

All of the above charges have to be handled with care because the integrals defining them do not 
generally converge for configurations with no n- vanishing degree. In jT2| a related field theory 
was studied and it was pointed out that the Noether charges can be related to moments of the 
vorticity. In our cases the moments 

J[^] = -vr / [xl + xD^ (2.22) 

and 

Pi[ilj] = 27i X2'y and P2M = -27r / X17 (2.23) 

are also conserved during time evolution according to ()1.11|) and can be obtained from the Noether 
charges J, Pi and P2 by integration by parts. For some configurations the integrals defining J 
and Pi are convergent even when those defining J and P are not. The Galilean symmetry also 
implies a conservation law, but it will not be required in this paper. 
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3 The n = 2 saddle point and its unstable mode 



Imposing extra symmetry is the easiest way of finding static solutions of the Ginzburg-Landau 
equations. The largest invariance group one can impose on configurations of degree n is a group 
containing reflections and combinations of spatial rotations R{x) by an angle x with suitable 
phase rotations. More precisely we define the subgroup 

Rn = CK {(i?(x),e*"^) G S0{2) X U{l)\x e [0,2n)} (3.1) 

of the symmetry group ()2.17|) . Configurations invariant under this group are of the form 

Mr,0) = fn{r)e'''', (3.2) 

where is real and satisfies the boundary condition 

/„(0) =0 for n 7^ and lim /„(r) = 1. (3.3) 

r— >oo 

The Ginzburg-Landau equation implies the following ordinary differential equation for /„ 

As proved for example in [T] , this equation has a unique solution satisfying the boundary condi- 
tions ()3.3|) for each n, and in the following we use fn to denote that solution. Near r = and 
for large r the equation can be solved approximately 

fni'f') ~ ^nf'^ for small r, 

fn{r) ~ 1 - TTT for larger, (3.5) 
2r^ 

where An are real constants. However, no exact expression for in terms of standard functions is 
known. For |n| = 1 the resulting vortex configuration is a stable solution of the Ginzburg-Landau 
equation. For |?t,| > 1 the solution is an unstable saddle point. Bounds on the number of unstable 
modes were derived by Ovchinnikov and Sigal in |2] . It follows from their analysis that for n = 2 
there is precisely one unstable mode. Since this mode plays an important role in our analysis, 
we describe and compute it explicitly. For further general background we refer the reader to j2j. 
The required unstable mode is an eigenvector of the Hessian of the Ginzburg-Landau energy 

RessE[^] = (^ JgV (3.6) 

For the saddle point ■?/'2 the Hessian takes the form 

„ 1 / -A + 2p2 - 1 P2e^'' \ 7^ 

Hessi?[^2] = 2( p^,~.e -A + 2^^ - 1 J ' ^^^^ 
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where p2 = |V'2p- In terms of the notation (|2.8|) the eigenvalue equation 

2Hess£;[V'2]f= (3.8) 

is thus equivalent to 

- + (2p2 - l)e + P2e'''^ = Ae, (3.9) 

which can also be obtained by linearising the Ginzburg-Landau equation around the saddle point 
'?/'2- If ^ solves ()3.9|1 then the leading term in the difference E[iIj2 + C,] — -E'['?/'2] is the second order 
term 



ImCO, (3.10) 



where we used the inner product 



{i,V) = l ji^V + vOd'x. (3.11) 

To study the eigenvalue problem ()3.9|1 and to compute the variation in the energy it is best 
to expand the 6'-dependence of ^ into a Fourier series. The term p2e^*^ in the equation leads to 
a coupling between Fourier modes e^^^ with k differing by 4, but all other modes decouple. As 
explained in the unique eigenfunction with a negative eigenvalue is of the form 

^{r,e) =u{r) + v{r)e^'^. (3.12) 

Inserting this expression into ()3.9p leads to coupled equations for the radial functions u and v 

1 d f du\ 

r— \ + [2p2 — l)u + = Xu 



r—] + [— + 2p2-l]v + p2U = Xv. (3.13) 



r dr \ dr J 
1 d / dv\ /16 
r dr \ dr J \ 

In terms of u and v the second variation of the Ginzburg-Landau energy for a fluctuation ^ of 
the form (j3.12p about the stationary point il)2 is 



oo 







rdr, (3.14) 



where we have assumed regularity of u and v at the origin and exponential decay at infinity in 
the integration by parts. This assumption will be justified below. Note that, in terms of the 
functions u and f , 

(f,f) = 27r/ {u^ + v'^)rdr. (3.15) 
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The eigenvalue (|3.13|) was first studied numerically by Dziarmaga in using a shooting 
method. Near r = the leading terms for a regular solution are 

u{r) ~ 1 — and v{r) ~ ar^, (3.16) 

where a is an undetermined parameter. For large r, it is convenient to write the equation in 
terms of the linear combinations = u + v and u_ = u — v: 



1 d f du 

r dr \ dr j \ / r 



Now the equations decouple for large r and one deduces the asymptotic form 

exp(-^5T2r) exp(v^^T2r) 



~ ^exp(-v^r)^^exp(v^r)_ ^^^^^^ 



Solving the eigenvalue problem with a shooting method means determining the two parameters 
a and A so that the solutions and m_ both decay exponentially for large r i.e. B = D = 0. 
Such a parameter search in a two-dimensional parameter space is tricky. We start our search by 
assuming that v is negligible compared to u and take the trial function ^o('^, 0) = uo{r) where 

Uo(r) = sech ( ] , (3.20) 



4 



One then finds that 



E'^'\{o) = -0.378 X ^(fcfo), (3.21) 

showing that A < —0.378 and hence that the value of —0.168 for the eigenvalue given in JH] is 
incorrect. In our search, starting with ,^0; we find the eigenvalue 

A = -0.41869, (3.22) 

with eigenf unctions u and v shown in figure 1. Note that u is qualitatively very similar to 
Uq. However, v is non- vanishing and negative. This is not surprising. Looking at the energy 
expression (j3.14j) we note that it is energetically favourable for the functions u and v to have 
opposite signs. 



4 Unstable manifold and truncated dynamics 

The unstable manifold of a saddle point is the union of paths of steepest descent from the saddle 
point. To define paths of steepest descent, one requires both a metric and a potential. Suppose 
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Figure 1: The unique negative eigenmode o f the Hessian at the saddle point il)2- The plot shows 
the radial functions in the parametrisation \3.1^) . 



are local coordinates on a finite-dimensional manifold with potential energy V and the metric 
is represented by the matrix Qij. Then the paths of steepest descent are solutions to the gradient 
flow equations 



dr 



dV 
dx^' 



(4.1) 



with x*(r) approaching the relevant saddle point as r ^ — oo. It is important to keep track of 
the metric even when it is fiat. For example, gradient flow on with the standard fiat metric 
takes the following form in complex coordinates z = Xi+ 1x2'- 



dz 
d^ 



— = -2 



dz ' 



(4.2) 



In the field theory we are considering here the potential energy is the Ginzburg-Landau 
potential energy functional (jl.2j) and the metric is the fiat metric. Since we work with complex- 
valued fields, tangent vectors are also complex-valued functions and the metric is the inner 
product of tangent vectors defined in ()3.1H) . The gradient flow equation is the following non- 
linear heat equation 



dip 



(4.3) 



Here r is an auxiliary "time" variable, which we distinguish notationally from the time variable 
t used in the Gross-Pitaevski equation. 

We have solved this equation numerically on a grid of size 100 x 100, corresponding physically 
to the square [—20, 20] x [—20, 20] in the XiX2-plane. To generate the gradient flow curve we 
start with the saddle point configuration ip2 and add a small perturbation dip = e^, where ^ is 



11 



the negative eigenmode (|3.12j) found in the previous section. At the boundary of the lattice we 
impose the Dirichlet boundary condition iplr, 9) = e^*^. The resuhing gradient flow curve does 
not depend significantly on e provided it is small enough. When e = 0, the discretisation effects 
provide sufficient perturbation to make sure that the gradient flow curve moves away from the 
saddle point ip2, but the initial flow is very slow. It is therefore numerically convenient to use 
a non-vanishing value of e. In figures 2 and 3 we show snapshots of the field and of the energy 
density during the early stage of the gradient flow, generated by starting with 6tp = 0.0005^. 
During the gradient flow, the unstable n = 2 vortex splits into two n = 1 vortices which drift 
apart. The saddle point solution ip2 has an energy density which has a degenerate maximum 
on a ring. The splitting process begins with the development of two local maxima, as shown in 
configuration A in figure 3. As the splitting process continues, these local maxima become two 
distinct vortices. Note, however, that during the early phase of the splitting process, the vortices 
overlap and deform each other significantly. This is clearly visible in configuration B in figure 
3. Only in the later stages of the splitting process do the vortices resemble two standard n = 1 
vortex solutions. 
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Figure 2: Vortex configurations in the XiX2-plane during the early stages of the gradient flow 
from the saddle point ip2 ■ 



In our numerical simulation of the gradient flow we generated 700 configurations. The most 
convenient way of labelling these configurations is in terms of their Ginzburg-Landau energy 
rather than by the configuration number or the gradient flow "time" r, which have no physical 
significance in the Gross-Pitaevski model. The finite lattice used in the computation provides a 
natural renormalisation of the energy of the vortex configurations, assigning the value i?rcn['^2] = 
Knax = 31.982 to the saddle point configuration '?/'2- The configuration A in figures 2 and 3 has 
the energy V = 31.955 and the configuration B has the energy V = 30.257. 
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Figure 3: Energy density in the XiX2-plane 
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the vortex configurations A and B of figure 2. 



There is another way of parametrising the gradient flow which we need to address, namely 
the parametrisation in terms of vortex separation. Vortex separation is used as a parameter 
for pinned vortex configurations in the work of Ovchinnikov and Sigal. In view of our later 
comparison with that work we therefore make a few comments on it here. The pictures of the 
vortex configurations and energy densities shown in figures 2 and 3 illustrate that the notion of a 
separation is not well-defined during the early stages of the gradient flow. The best one can do is 
to define a functional of vortex configurations which reproduces the distance between the maxima 
of the energy density (or, equivalently, between the zeros of the field ip) for well-separated vortices 
and which is computationally convenient. In the following we use the vorticity ()2.(jj) instead of 
the energy density. The vorticity looks qualitatively very similar to the energy density during all 
stages of the gradient flow and in particular it is strongly peaked around the zeros of the field if) 
for well-separated vortices. The advantage of using the vorticity rather than the energy density 
is that the total vorticity is conserved whereas the total energy changes during the gradient flow. 
To motivate our approach further, note that the configurations generated during the gradient 
flow break the invariance group R2 of the saddle point to the group Z2 x Z2 generated by 
refiections about the xi and X2 axes, each combined with the complex conjugation of ip. As a 
result the vorticity 7 of all conflgurations is reflection symmetric about two orthogonal axes. For 
our gradient flow the reflection axes are the coordinate axes, and the vortices separate along the 
X2-axis because we have broken the rotational symmetry R2 of the saddle point solution explicitly 
by using the perturbation ^ ()3.12j) . A different perturbation would have led to an isomorphic 
reflection symmetry about a different set of orthogonal axes, and a different separation direction. 
We use the reflection axis orthogonal to the separation direction to divide M? into two half spaces 
Hr = {{xi,X2) G M'^\x2 > 0} and Hi = {{xi,X2) G M^|a;2 < 0}. Then we deflne the separation 
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functional as 

a [■?/'] = / X2'jd'^x— / X2'jd'^x = 2 / X2'jd'^x. (4-4) 

As expected, the separation functional ()4.4|) increases monotonically during the gradient flow, 
but one interesting feature is that it assigns a non-vanishing value to the "coincident" vortex 
configuration i\}2- This is a familiar feature of separation parameters used for other solitons, such 
as monopoles in non-abelian gauge theory jTHj- The point is that solitons lose their separate 
identity when they overlap and separation ceases to be a meaningful concept. In particular 
there is no reason to insist that a separation parameter vanish for a particular configuration. 
Numerically, we find the minimal separation to be a^:,^ = a[^/'2] = 4.815, which corresponds 
roughly to the diameter of the ring on which the energy density of 1^2 is maximal. For the 
configuration A in figures 2 and 3 we find a = 4.896 and for configuration B we find a = 6.406. 
The lowest energy configuration we generated has the energy V = 24.405 and separation a = 
15.570. 

For the construction of the unstable manifold we return to the parametrisation of the gradient 
flow curve in terms of renormalised energy. Thus we obtain a family of fields "^{V; x), labelled by 
the value V of the renormalised Ginzburg-Landau energy. According to the general prescription 
of JU] we should now act with the full symmetry group ()2.17|) to generate a family of saddle 
points and a family of gradient flow curves, and use their union as the collective coordinate 
manifold for the truncated dynamics. However, we are only interested in the relative motion of 
the two vortices in the fission process and not in their centre-of-mass motion. The latter can be 
found by applying Galilean boosts ()2.18|) to the entire configuration. Furthermore, we exclude 
collective coordinates which change our boundary condition. Both phase rotations and spatial 
rotations change the field "at infinity" (or at the boundary of our lattice), but the subgroup 
i?2 defined in p.lj) respects the boundary condition. It also leaves the saddle point solution ip2 
invariant but it acts non-trivially on the gradient curve emanating from il)2- Thus we obtain a 
family of fields 

^{V,x;r,9) = ^{V;r,e-x)e'''' (4.5) 

depending on two collective coordinates V, x- The former measures the energy and the second 
labels the spatial orientation of the two-vortex configurations generated during the gradient 
flow. Since the saddle point configuration with energy V|„ax does not depend on the angle x, the 
collective coordinate manifold is topologically a plane, with (Knax — V, x) as polar coordinates. 
We denote this manifold by Ai2- The range of the angle x is [0, tt) since the configurations 
"^{V;x) generated during the gradient flow are invariant under a rotation by tt (which is the 
product of the reflections at the coordinate axes discussed above). Note that A42 C €2- 

It is now a simple matter to compute the restriction of the Gross-Pitaevski Lagrangian 1)2.11) 
to 7^2- We allow the collective coordinates V and x to depend on time and insert the fields 



14 



'^{V{t),x{t);r,9) into (EH). The result is 

L[^{V, x; ■)] = -\imV) - J{V))x - Ig{V)V - V, (4.6) 

where we have written Q{V) and J{V) for the Noether charges ()2.19p and ()2.20|) evaluated on 
^{V,x;x), i.e. 

QiV) = J \^iV,x;x)\'d'x (4.7) 

and 

J{V) = j \m(^{V,x;x)^{V,x;x)^ £x. (4.8) 

The integrals defining Q and V are both divergent, but the combination 

K{V) = J{V) - 2Q{V) (4.9) 

occuring in the Lagrangian is finite. To see this note that for the saddle point configuration iIj2 
one computes J = 2Q, so that K{V^^^) = 0. As the configuration evolves away from the saddle 
point K grows but it remains finite for any finite value of V. This follows from the continuity 
of the semigroup defined by the evolution equation 14. HI in a suitable norm, see [T^ and 
and from the continuity of the funtionals J and Q ()2.19|) and ()2.2()j) with respect to that norm. 
Finally, the function G{V) is 



GiV) = J lm(^^iV,x;x)^{V,x;x)^ £x. 



(4.10) 



Neither K nor G depends on x because of the rotationally invariant integration measure and 
because the integrands are independent of the phase of \E'. The term G{V)V is a total time 
derivative which does not affect the equations of motion. We therefore omit it in our final 
expression for the induced Lagrangian L2 on A^2: 

L2 = ^K{V)x-V. (4.11) 

The Euler-Lagrange equations are very simple because our truncated system, like the original 
Gross-Pitaevski dynamics, is rotationally invariant and has a conserved energy. Thus V remains 
constant during the time evolution, and the angle x changes according to 

xiV) = ^ (4.12) 

It follows from the constancy of V that x is also constant during the time evolution. We introduce 
the abbreviation 

uiV) = xiV) (4.13) 
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for the rate of change of the angle x- The dependence of this angular velocity on V is the main 
dynamical information we extract from our truncated dynamics. It is is given by the differential 
quotient 

-(n = 2§. (4.14) 

Numerically, we compute uj by evaluating finite difference quotients AV^/A/T of successive con- 
figurations. Our results are displayed in figure 4. In the early phase of the gradient flow (roughly 
the first 100 configurations) both V and K change very slowly and accurate numerical compu- 
tation of the difference quotient is difficult. We have omitted these difference quotients from 
our plot. Since the energy V hardly changes during this part of the gradient flow, this omission 
makes no visible to difference to the graph ojiV) shown in figure 4. The next section is devoted 
to a detailed discussion and interpretation of our results. 

0.25 
0.2 
0.15 

to 

0.1 
0.05 



'2A 25 26 27 28 29 30 31 32 

V 

Figure 4: The solid line shows the angular velocity uj of two-vortex configurations as a function 
of the energy V. The dependence of uj on the energy predicted by the Kirchhoff-Onsager law is 
plotted with a dashed line. The value of uj„,^^ given in ( 1.5. .9)) is marked with a "x 




5 Orbiting vortex pairs 

Our two-dimensional Lagrangian dynamical system {M.21 -^2) is designed to approximate the slow 
interactive dynamics of two vortices moving according to the Gross-Pitaevski equation ()1.11|) . 
The model predicts that the two vortices will orbit each other, with the angular frequency 
depending on their total energy as shown in our figure 4. While the qualitative behaviour of 
the two vortex system is not surprising and familiar from vortices in ordinary fluids and also 
vortices in the gauged Ginzburg-Landau model jTH], the precise prediction of the variation of the 
angular frequency with the energy for both overlapping and well-separated vortices is new. In 
Dziarmaga computed the rotation frequency of an overlapping vortex pair by linearising the 
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Gross-Pitaevski equation around the saddle point solution ijj. The linearised solution is 

-^{t, x) = + ee-'^\u{r) + e^*%(r)), (5.1) 

where A and {u,v) are the eigenvalue and eigenfunctions of the eigenvalue problem (j3.13|) . The 
interpretation of the linearised solution in terms of rotating vortices requires some care. As we 
saw in section 3, the radial eigenfunction v is much smaller than u. Thus, as a first approximation, 
we may think of (jS.ip as a field whose magnitude depends on r according to eu{r) and whose 
direction changes with angular frequency A, all superimposed on the static n = 2 solution. The 
effect of such a superposition is a configuration with two zeros (whose separation depends on 
e) which rotate around one another with angular velocity = —A/2. The factor 1/2 arises 
because for ip2 a phase rotation by a is equivalent to a spatial rotation by a/2. With our value 
for A, we thus find = -A/2 = 0.20935. 

In the limit V — > V^^^ our numerically computed function uj{V) approaches a value which is 
close to but not exactly equal to Un^. The discrepancy arises because we neglected the radial 
eigenfunction v in our interpretation of the linearised solution. We can establish a precise relation 
between the eigenvalue A and the limiting angular frequency 

a;,^,, := lim uiV) (5.2) 

y ' vmax 

as follows. Using the definition (j4.9|) and noting that both V and K are originally defined along 
the gradient flow curve parametrised in terms of the auxiliary variable r fj4.3p we write 

dV dV ,dK , , 

Recall that we computed our gradient flow curve by starting at r = with the configuration 
V'2 + e^, where ^ is the eigenmode ()3.12p and e is small (in practice e = 0.0005). At that point 
the gradient of the Ginzburg-Landau energy functional is 

= ^Aee, (5.4) 

where we have used that the first functional derivative of E vanishes at V'2 and that ^ is an 
eigenfunction of the Hessian, i.e. it satisfies ()3.8p . Hence from ()4.3p 



and thus 

dV 
dr 



T = 



drip\,^o ~ -Ae^ and 5^^|,^o ^ -Ae^ (5.5) 

6E 6E 

^[V^2 + e^] dripi^^o + ^[^2 + e^] dripi^^Q d^x 

-e^X' J Ud^x 

-e^A^ ■ 27r / {i? + v"^) rdr, (5.6) 
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where we used (I5.5p and inserted the expression (j3.12p for ^. 

It is easy to check that any i?2-invariant configuration ()3.H) and hence in particular 11)2 is a 
stationary point of the functional K[il)\ fl4.9p . Thus |^(V'2) = and by a similar calculation to 
the one given above we have 



dK 
dr 



T=0 



/S K S K 



J Sip Sip Sip Sip 
= -2e^\ J {Im^de^ - 2^^ d'^x 

= Ae'^X ■ 271 j {u^ - f 2) rrfr, (5.7) 

where we used that S'^KjSi\P' = S'^K/Sip'^ = 0. The approximate equalities become exact in the 
limit e — > 0. Thus we find the following exact expression for the angular frequency w,^^^ ()5.2p in 
terms of the eigenvalue A: 

2 J (m^ — v^jrdr 

In view of the numerical difficulties in computing uj near the saddle point this formula is very 
useful in practice. Evaluating it numerically we find 

tu^,, = X 1-0200 = 0.21354. (5.9) 

This value provides a check on our computation of ioiy^ near the saddle point and is marked 
with a "x" in figures 4 and 5. 

For well-separated vortices we can compare our results with those obtained by Neu in O E] 
via scaling techniques or those obtained by Ovchinnikov and Sigal in [S^ based on pinned vortex 
configurations. Both find that well-separated vortices are governed by the Kirchhoff-Onsager 
law. According to that law two vortices separated by a distance a have an interaction energy 
given by 

y(a) = -27rlna + C, (5.10) 

where C is an arbitrary normalisation constant, and orbit each other with angular frequency 

4 

t^Ko = (5.11) 
It follows that the rotation frequency is given as a function of the energy via 

a7Ko(^) =C'e^, (5.12) 
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where C = 4exp(— C/vr). The constant C (and hence C) is determined once the dependence 
of the energy on the vortex separation is known. To make contact with the Kirchhoff-Onsager 
relations we thus have to make use of the separation functional introduced for our vortex con- 
figurations in (|4.4j) . We pick a vortex configuration consisting of clearly separated vortices and 
compute its energy to be ^ = 27.253 and its separation to be a = 10.012. Inserting these values 
into ()5.10p fixes the value of C and hence C. We have plotted the resulting function ()5.12|) 
in figure 4. It is important to be clear about the interpretation of the two curves in figure 4. 
For well-separated vortices we can meaningfully associate a separation parameter to a vortex 
configuration. The Kirchhoff-Onsager relation gives the angular frequency of point-vortices with 
that separation. The fact that for low energy configurations consisting of well-separated vortices 
the Kirchhoff-Onsager relation agrees with our results shows that well-separated vortices may 
be treated as point- vortices, in agreement with the results of Ovchinnikov and Sigal in |[H|. The 
two curves in figure 5 differ for configurations near the saddle point, not so much because the 
predictions of our model disagree with those of the point vortex model but because the two 
models are no longer compatible in that regime. 

From the point of view of the Kirchhoff-Onsager relation, the dependence (IS.llj) of the angular 
velocity on separation is more fundamental than the dependence on the energy. In figure 5 we 
plot the angular velocity as a function of the separation parameter ()4.4j) and compare it with 
the simple Kirchhoff-Onsager prediction. Again the two plots agree for large separations, and 
the comments made above on the comparison near the saddle point apply here, too. However, 
it is worth emphasising that the Kirchhoff-Onsager relation predicts a divergence of the angular 
velocity when the separation tends to zero. Our results show that for the Ginzburg-Landau 
vortices this divergence is removed, essentially by cutting off small separation parameters. This 
is yet another example of a soliton model "regularising" the singularities of a point-particle 
model. 

6 Discussion and Conclusion 

The unstable manifold method used in this paper to approximate the dynamics of interacting 
vortices in the Gross-Pitaevski model allows one to analyse the relative motion of two vortices 
at arbitrary separation. When the vortices are overlapping our results reproduce those found 
with linearisation methods and when the vortices are well-separated our results agree with the 
Kirchhoff-Onsager law. Our method gives a unified description of both regimes, and successfully 
interpolates between them. 

When vortices evolve according to the Gross-Pitaevski equation they will in general also emit 
radiation. This effect is not captured in our approximation. However, it is possible to compute 
radiative corrections to the vortex motion predicted by our approximation using the methods 
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Figure 5: The solid line shows the angular velocity uo of two-vortex configurations as a function 
of the separation a. The dependence of uo on the separation predicted by the Kirchhoff-Onsager 
law is plotted with a dashed line. The value of 00^^.^. given in \5.y\) is marked with a "x 

developed by Ovchinnikov and Sigal in In particular, Ovchinnikov and Sigal studied the 
radiative corrections to the orbiting motion of a vortex pair and found that, for large vortex 
separation a, the power emitted in the radiation is proportional to a~^. As a result of this energy 
loss, the orbiting vortices drift apart slowly as they orbit each other. However, since the power 
loss is small, the change of the separation during one period of rotation is also small (much less 
than 1 %). Combining these results with our calculations, we arrive at the following qualitative 
picture of vortex fission in the Gross-Pitaevski model. Suppose the n = 2 saddle point solution 
is perturbed and allowed to evolve according to the Gross-Pitaevski equation. The perturbed 
vortex configuration rotates and emits radiation. The radiation carries away energy and the 
vortex configuration adjusts by finding a configuration of lower energy. The most efficient way 
of doing so is by following a path of steepest (energy) descent. As we have seen, the n = 2 
vortex breaks up into two n = 1 vortices along the path of steepest descent. Thus we can use 
our truncated dynamical system {M.21 -^2) and expect the vortex pair to orbit with the angular 
frequency depending on the energy according to ()4.12j) . The continued rotation results in further 
emission of radiation and loss of energy, to which the vortex pair adjusts by drifting further apart. 
The fission process is thus governed by the combination of two effects: the radiation carries away 
energy and the vortex configuration adjusts by sliding down the unstable manifold of the n = 2 
saddle point. It would be interesting to make this picture more precise by studying the radiation 
effects quantitatively. This should probably be done in conjunction with a more careful study of 
the conservation laws in the theory, similar to the investigation for gauged vortices in As 
mentioned briefly in section 2, the conservation of angular momentum can also be expressed as 
the conservation of the moment ()2.22|) of the vorticity. Just looking at vortex motion during the 
fission process it seems that the moment ()2.22|) increases as the vortices drift apart. It would 
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be interesting to understand how the radiation makes up for this change, possibly by carrying 
negative vorticity off to infinity. 

We would like to emphasise that the work described here is the first application of the 
unstable manifold method to soliton dynamics in more than one spatial dimension. As such, it 
contains a number of useful general lessons. At first sight it may not seem sensible to use one 
non-linear partial differential equation (the gradient flow equation) to approximate another non- 
linear partial differential equation (the Gross-Pitaevski equation). However, generally speaking, 
gradient flow is easier to understand and to compute numerically than Hamiltonian flow like 
that defined by the energy- conserving Gross-Pitaevski equation. One reason for this is that the 
typical time scales of the two evolution processes are very different. Our discussion of the vortex 
fission process illustrates this point. The solution of the gradient flow equation over a relatively 
short CPU time (roughly 700 steps) maps out all the interesting configurations between the n = 2 
saddle point and the well-separated vortices. The fission process according to the Gross-Pitaevski 
equation, by contrast, is expected to take several orders of magnitude longer, making it difficult to 
maintain numerical accuracy. The gradient flow equation provides an efficient way of mapping out 
those configurations which are relevant for slow (or low-energy) dynamics. At a more practical 
level, we found that the explicit study of the linear negative modes of the saddle point was 
essential for an efficient computation of the unstable manifold. We expect that this will become 
even more important for saddle points with more than one negative mode. Thus one could repeat 
the calculations performed here for the saddle point solutions ipn with vortex number n> 2. As 
explained in there is now more than one unstable mode. As a first step in constructing the 
unstable manifold of these saddle points one needs to compute the negative eigenmodes of the 
linearised Ginzburg-Landau equation explicitly. While the explicit construction of the unstable 
manifolds and the computation of the induced dynamics may be difficult, it would be interesting 
to understand general geometrical features of the unstable manifolds and use them to derive 
qualitative features of multi-vortex dynamics in the Gross-Pitaevski model. 

It remains an open problem to justify the unstable manifold approximation used in this paper 
analytically. The fact that our results agree with those obtained via more familiar approximations 
in limiting regimes provides an encouraging check. However, this agreement also raises a question. 
Near the saddle point our method agrees with results obtained by linearisation, and the small 
parameter controlling the validity of the approximation is the amplitude e in ()5.H) . As explained, 
this is related to the separation of the two zeros of the field tp. Far away from the saddle 
point, our approximation agrees with the point-vortex model implicit in the Kirchhoff-Onsager 
law. The small parameter controlling the validity of the approximation in this regime is the 
inverse separation of the zeros of if) (i.e. the inverse distance between the vortices). Since 
our method interpolates between the two approximations it is not clear which, if any, small 
parameter controls the approximation in the intermediate regime. More generally, proof of the 
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validity of the approximation, possibly along the lines of the proof of the geodesic approximation 
for gauged vortex dynamics in would be very desirable. A useful starting point for such 
an investigation may be the following natural relation between the gradient flow equation ()4.3p 
and the Gross-Pit aevski equation (jl.llj) . Both define flows in the space C (|2.4j) in terms of the 
Ginzburg-Landau energy. However, the gradient flow uses the inner product (j3.1H) whereas the 
Gross-Pitaevski equation is a Hamiltonian flow using the symplectic structure ()2.9|) . The inner 
product ()3.1H) and the symplectic form ()2.9|) are equal to, respectively, the real and imaginary 
part of the sesquilinear form 



One interesting manifestation of this relationship is that gradient flow trajectories and solutions 
of the Gross-Pitaevski equation are orthogonal with respect to the inner product ()3.1H) whenever 
they intersect. 
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